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[nr] ■ We investigate the propagation of accretion-powered jets in various types of 

■ massive stars such as Wolf-Rayet stars, hght Population III (Pop III) stars, and 

massive Pop III stars, all of which are the progenitor candidates of Gamma-Ray 
O . Bursts (GRBs). We perform two dimensional axisymmetric simulations of rela- 

^ ■ tivistic hydrodynamics taking into account both the envelope collapse and the jet 

propagation (i.e., the negative feedback of the jet on the accretion). Based on our 
hydrodynamic simulations, we show for the first time that the accretion-powered 
'i> . jet can potentially break out relativistically from the outer layers of Pop III pro- 

. genitors. In our simulations, the accretion rate is estimated by the mass flux 

^ : going through the inner boundary, and the jet is injected with a fixed accretion- 

to-jet conversion efficiency 77. By varying the efficiency rj and opening angle 6 op 
O ' for more than 40 models, we find that the jet can make a relativistic breakout 

from all types of progenitors for GRBs if a simple condition 77 > 10~^{9op/8°)^ is 
satisfied, which is consistent with analytical estimates. Otherwise no explosion 
or some failed spherical explosions occur. 
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Introduction 



The link between nearby Gamma Ray Bursts (GRBs) and peculiar Type Ib/c super- 
novae (or hypernovae) ambiguously shows that some populations of GRBs are born from 
the catastrophic death of massive stars. Observations of host galaxies of GRBs also lead to 
the general cons ensus that GRBs are generated preferentially in low metallicity star-forming 
regions (see e.g. iModjaz et al.ll2008l ). The stellar evolutions at low metallicity are believed 
to suppress the mass loss, so that the star can maintain its own angular momentum. As 
a result, the iron core of the se stars would be rapidly spinning (see e.g. lYoon fc Langer 
20051 : IWoosley fc Hegerl 120061) . which would be a necessary c ondition for producing GRBs 
in the coUapsar ( WooslevI 19931: MacFadven fc Woosley 1999 ) or magnetar models (see e.g. 



Thompson et al. 



2004 



Metzger et al.ll201ll ). According to these facts, it is widely recog- 
nized that rapidly rotating Wolf-Rayet stars in low metallicity regions are the most favored 
progenitors for GRBs. 

The first stars (hereafter Pop III) also potentially cre ate GRBs. The first stars are 
supp osed to be formed with a huge mass (M > lOOMj 







2OO2I ) and a rapid rotation with nearly breakup speed (IStacv et al.l I2OIII). The gravita- 



flAbel et al. 



tional collapse of these stars would result in a black hole formation (INakazato et al.l 12006 



Suwa et al. 



2002; Bromm et al. 



2007, 



2009 



Sekiguchi fc Shibata]|201ll) and potentially the c e ntral engine activit y 



faeger et al.lboosi iMeszaros &: Reeslboioi iKomissarov &: Barkwiboioi Isuwa fc lokalboilh . 



If the primordial gas is ionized by radiation from first-generation me tal-free (Pop III.l) s tars 



the subsequent metal-free (Po p III. 2) stars would b e less massive ( iBromm et al.l 120091 ) and 



outnumber the Pop III.l stars fide Souza et al.ll201ll ). It has also been recently discussed that 



the radiative feedback could reduce the mass o f the first star via the HII region breakou t 



and the photoevaporation of the accretion disk f lMcKee fc Tanll2008l : iHosokawa et al.ll201ll ). 
Although the reduced mass < lOOM© is significantly lower than previously thought, the light 
Pop III stars still have large enough mantles for the formation of a black hole. Thus, even in 
these cases, the central engine could operate as a result of the core collapse in the standard 
coUapsar model. The Pop III GRBs and their afterglows are detectable in principle up to 



(Lam 


3 & Reichart 


2000 


„| . . , , ^,^w.- 

Ciardi & Loeb 


2000; 


loka 


2003 


2005; 


Inoue et al. 


I2OO7 




Toma et al. 


201] 





Iokall2003l : IGou et al.l 120041 : lloka fc Meszaros 



However, even if the central engine successfully operates, it is still a matter of debate 
whether the jet can produce GRBs or not. One of the main obstacles for producing GRBs is 
the stellar envelope that may prevent the jet propagation. If the central engine turns off well 
before the jet head reaches the stellar surface, all of the jet matter undergoes dissipation by 
the reverse shock wave and it will eventually expand spherically. In addition, the outfiow 



- 3 - 



is contaminated by a huge amount of baryons, so it is naturally expected that its velocity 
becomes non-relativistic and never create a GRB. With this expectation, iMatznerl (j2003[ ) 
constrained the progenitors of GRBs assuming that the lifetime of the central engine is 
comparable to the observed duration of the prompt phase of GRBs. He concluded that only 
compact carbon-oxygen Wolf-Rayet stars satisfy the condition for producing GRBs, while 
very massive stars such as Pop III stars are not suitable. 



On the contrary, ISuwa fc lokal (120 111 ) recently pointed out that the jet breakout is 
possible even if the Pop III star has a supergiant hydrogen envelope without mass loss, 
thanks to the long-lived powerful accretion of the envelope itself. They analytically showed 
that the jet successfully penetrates the Pop III as well as compact Wolf-Rayet stars if the 
envelope continues to fall in a black hole and the accretion-to-jet conversion efficiency is 
larger than a certain level. 

However, it is not trivial to determine whether the envelope can continue to fall in and 
accrete onto black holes or not. Generally, the core collapse produces rarefaction waves, 
whi ch propagate outwards through the envelope and induce the infall of the stellar envelope 



see 



Nagakura et al.ll201ll ). But some portions of the envelope cease to fall due to the jet 



propagation when the central engine begins to operate. Although almost all matter could 
accrete from the equatorial regions, the feedback would affect the accretion rate if the jet 
opening angle is la rge. The jet feedback to the accretion has not been take i i into account 
in pr e yious studies 



2010; Suwa fc loka 



lMatznerll2003l : Ijaniuk fc Progall2008l : Kumar et al.ll2008al : iLindner et al. 



201ll ). Since these processes are supposed to be complex and strongly 



non-linear phenomena, hydrodynamic simulations with both accretions and jet propgations 
are strongly required. 



On the other hand, a larg 



mant l e have been carried out ( MacFadven fc Woos. 
200?! 



Tominaga et a. 



2007 



e number of numerical studies on jet propagations in the stellar 



Lazzati et al. 



20091: 



19991: 



Mizuta et a" 



Mizuta fc Aloy 



2009 



2006[lMorsonv et al 



Mizuta et al 



2011 



Nagakura et al.ll201ll ). However, almost all works assume that the jet is injected with a 
constant energy flux from a certain radius of the inner boundary. Although we do not know 
the mechanism of the central engine, it is naturally exp ected that the jet lumi r iosity would 



corre l ate with the accretion rate in one way or anoth er (jPi Matteo et al 



2002[ IProga et al. 



20031 : lMcKinnevll2006l : IZalamea fc Beloborodovll201l[ ) (but see also flMetzger et al.ll201lh for 
a magnetar model that does not depend on accretion). We also wonder whether the jet 
production could be ceased by the reductio n of the accretion d ue to the negative feedback 
as mentioned above. In the previous study, iMacFadyen et al.1 (120011 ) demonstrated the jet 
propagation with the jet luminosity as a function of the accretion rates. However, they cal- 
culated the jet propagation and the fall back process separately, and can not address the jet 
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feedback process adequately. In another previous study, iMorsony et al.l (120 lOf ) investigated 
time variable jet injections, but the luminosity and time variability are determined by hand. 
In order to investigate the jet propagation feedback on the accretion process, it is necessary 
to perform numerical simulations in both the collapsing phase and the jet propagation phase 
at once. 

Motivated by these facts, we perform two-dimensional axisymmetric hydrodynamic sim- 
ulations for the envelope collapse and the jet propagation in a single computation. The 
purpose of this study is to clarify whether the forward shock wave successfully propagates 
and breaks out from the various types of the stellar progenitors for GRBs or not, taking into 
account the jet feedback process. We survey the parameter space of the jet opening angle 
and the accretion-to-jet conversion efficiency, and discuss how these key quantities affect the 
jet dynamics to obtain the simple analytical criteria for the GRB production. This paper 
is organized as follows. In Section 2, we describe the models and methods in this paper. 
Then, our results will be presented with detailed analyses in Section 3. Finally, we discuss 
our findings and conclude the paper in Section 4. 



2. Numerical Methods and Models 

The numer i cal co des employed in this paper are essentially the same as those used in 



Nagakura et al.l ( 1201 ll ). in which all the details about our numerical codes and various test 



calculations are presented. Here we briefiy summarize the methods and setups in this study. 

Our numerical code solves the relativistic hydrodynamic equations with a weak gravi- 
tational field. The self gravity is included in the weak field approximation of the Einstein 
equation. It should be noted that, due to our computational limitations, we cut the inner 
portions of the star from a certain radius. The gravity in this region is added as that of a point 
mass at the center by integrating the mass fiux at the inner boundary. Based on the above 
assumptions, the gravity is solved by using MICCG methods. The hydrodynamical parts are 
solved by using the so-called central scheme, which guarantees good ac curacy even if the fiow 



involves strong shock waves and the fiow velocity is highly relativistic (IKurganov &: Tadmor 



20001 : iNagakura fc Yamadall2008l ). We use the PPM interpolation method and TVD Runge- 
Kutta time integration which achieve second order accuracy in both space and time. In 
this study, we adopt the 7-law equation of state (EOS), p = (7 — l)poe with p, 7 = 4/3, 
Pq and e being the pressure, adiabatic index, rest mass density and specific internal energy, 
respectively, in all our computations. 



In this paper, we adopt three representative stellar progenitors, which are (1) Wolf-Rayet 
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star (16TI in flWooslev fc Heeerl 120061 ) . hereafter WR), (2) light Pop III star AOMq, which 



is the metal free pre-supernova model calculated by IWoos 



and (3) massive Pop III star 915Mq ( lOhkubo et al.ll2009l ) (hereafter mpop3). The density 



ey et al.l (|2002[ ) (hereafter IpopS), 



profile of each stellar model is displayed in Figure [H and the stellar mass and radius are 
also summarized in Table [H Since this study is the first attempt for the accretion-powered 
jet propagations, we neglect the stellar rotation for simplicity although these progenitors are 
supposed to spin rapidly to operate the centr al engine. Note that we use the same approach 
as in the previous study (jSuwa fc Iokall201ll ). that the accretion-to-jet conversion efficiency 
parameter absorbs this uncertainty. More detailed studies for the effects of rotation will be 
presented in the forthcoming paper (INagakura et al.ll2012l ). 



We map the spherical symmetric progenitors into two-dimensional grids in spherical 
coordinates. The computational domain for each stellar progenitor covers from a certain 
inner radius to the stellar surface. Although the inner boundary should be located in the 
vicinity of a black hole (around 10^~^cm), this is computationally very expensive, so that we 
set the inner boundary far from a black hole. For our reference models, the inner boundary 
for each progenitor is located at Rin ~ 10~^ x Rgtar, where Rgtar denotes the stellar radius 
(see Table [2]). According to this limitation, our discussions in the present paper are at 
the qualitative level. The outer boundary for each model is set in slightly outside from 
each progenitor. It locates at 4 x lO^^cm, 1.7 x lO^^cm and lO^^cm for WR, lpop3 and 
mpopS, respectively. For all models, the number of standard radial grid points is 500. Th e 
grid width is non-uniform and increasing in geometric progression ( iNagakura et al.ll201ll ). 
The innermost grid width is set as Ar = Rin/lO where R^ is the radius of the inner 
boundary from the center. Then, the rate of geometrical increase is determined so as to 
cover all computational regions with 500 meshes. The angular grid covers a quadrant of 
the meridian section (where we assume equatorial symmetry) and is uniform with 60 grid 
points. We employ an adaptive mesh refinement (AMR) technique in ord e r to d ecrease the 
computational cost. We deploy two levels of meshes as in INagakura et al.l ( 1201 ll ). where the 
resolution of the second level is 3 times finer in each direction than the first standard mesh. 
Thus, the angular resolution in AMR region corresponds to 0.5° for all models. The smallest 
radial grid width for WR model is 1.6 x lO^cm, while one for lpop3 and mpop3 are 3.4 x lO^cm 
and 3.4 X lO^cm, respectively. In order to check the dependence on the resolution, we also 
carry out finer AMR calculations. We check calculations with 5 times and 7 times finer 
AMR meshes for WR models (WRreso5 and WRreso7), and 7 times for other progenitor 
models (lpop3reso7 and mpop3reso7). For 5 times finer AMR meshes, the angular resolution 
corresponds to 0.3°. The smallest radial grid widths are 9.6 x lO^cm (WR), 2.0 x lO^cm 
(lpop3) and 2.0 x lO^cm, respectively. For 7 times finer AMR meshes, the angular resolution 
is (3/14)°. The smallest radial grid widths are 6.9 x lO^cm (WR), 1.5 x lO^cm (lpop3) 
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and 1.5 X lO^cm, respectively. As we shall see in Section 13. 6[ the numerical resolution is 
important to prevent baryon pollution by numerical diffusion (although it does not affect the 
macroscopic jet dynamics. See Section IXB]) . Note that in Table [31 the "AMR level" denotes 
the multiplying factor of the fine meshes. 

We assume that the central engine successfully operates and the well confined outflows 
are produced in the vicinity of a black hole. We inject the plasma in the radial direction 
through the inner boundary with an opening angle of several degrees. In this paper, we also 
assume that the jet luminosity depends only on the accretion rate which is estimated by the 
mass flows across through the inner boundary; 

M = -2n po{rin,ey{rin,9)un^smede, (1) 
Jo 

where and v"^ denote the location of the inner boundary and the radial velocity of flows. 
We inject an outflow with a luminosity, 

Lj,t = r]Mc^- (2) 

where rj and c denote the accretion-to-jet conversion efficiency parameter and the speed of 
light, respectively. The maximum conversion efficiency as a consequence of the accretion 
process is 5.7% for a Schw arzschild black hole and 42% for an extreme Kerr black hole 



f lShapiro fc Teukolskyi Il983l ). Thus, we choose the value of rj less than that. The conver- 
sion efficiency, jet opening angle, specific internal energy and radial velocity are varied in 
each model. Once these parameters are fixed, the density and pressure of injected jets are 
determined by using the relation, 

L,et = PoTv'ihT - lyAS (3) 

where h{= 1 + e/c^ +p/(Poc^)) and AS denote the specific enthalpy and the area of the 
injection surface, respectively. 



The collapse of the massive envelope is induced in the same way as in iNagakura et al. 



(120 111 ). In reality, the stellar envelope begins to fall after the arrival of a rarefaction wave 
that is generated by the inner core collapse. We mimic this situation by putting the radial 
gradient of all quantities to zero at the inner boundary except for the jet injected regions. 
At the beginning of the simulation, the break of the force balance at the inner boundary 
induces the infall of matter. Subsequently a rarefaction wave propagates outwards, inducing 
the infall as it reaches each point. It should be noted, however, that the stellar progenitors 
used in this p aper, especially t h e WR progenitor, are not exactly in dynamical equilibrium. 



in contrast to INagakura et al.l (1201 ll ). Accordingly, the outer parts of the envelope begin 
to move and the several artificial waves are observed during the simulations. We find that 
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they induce artificial explosions in some models (but only for non- successful shock breakout 
models). This is because the W R stellar model, for example, originally involves rotation 



[see e.g. 



Woosley fc Hegerll2006[ ). so that the centrifugal force works to sustain the stellar 
configurations. Since we artificially remove rotation, the envelope tends to infall and induce 
artificial compressions and bounces. We n ote that even th e origi nal 16TI model is not in 



exact dynamical equilibrium. Contrary to iNagakura et al.l (120111 ). we do not take special 



treatments for the initial stellar configurations here. 

We also investigate the dependence on the timing of the jet injection, since we still 
have few constraint on the starting time of the central engine. If the operation of the 
central engine is sufficiently late, the density profile of the stellar envelope is changed by the 
accretion. Thus, it is expected that the jet dynamics also depends on the timing of the jet 
injection. In addition, we would like to investigate whether the later jet can really accomplish 
the shock breakout since the large mass accretion may prevent the jet propagations in this 
case. Motivated by these facts, we initially let the stellar envelope spherically collapse, and 
then inject a relativistic jet. We carry out these simulations only in WR models because the 
enclosed mass at the inner boundary for our reference model is M;„ ~ 2M(7^, w hich may be 



lower than the critical mass of the black hole formation (jPemorest et al.l 120101 ). As a result, 
it is quite likely that the central engine does not operate for a while. On the other hand, 
the Pop III models have so much mass enclosed at the inner boundary (see Table [2]) that 
the central engine would begin to work soon after the collapse in our models. We prepare 
two models, WRM3 and WRM6, which inject the jet when the enclosed mass at the inner 
boundary reaches Mj„ = SM© and 6Mq , respectively. The corresponding retarded time of 
injection for each model is tiate = 7.47s and 26.93s, respectively. The radial density profiles 
just before the jet injection are displayed in Fig [2] for these models. 

All of our models used in this paper are summarized in Table [21 We prepare a reference 
model for each progenitor with the following parameters. The inner boundary is located at 
Rin ~ 10~^ X Rstari where Rstar denotes the stellar radius. The accretion-to-jet conversion 
efficiency (77), half opening angle of outflows (6'op), injected Lorentz factor (Fj), injected spe- 
cific internal energy (e^) are set as 77 = 10"'^, 6op = 9°, Tj = 400 and ej = 10~^, respectively. 
Note that these models assume that the injected jet has already reached to the terminal 
Lorentz factor. It may not be true for WR models since the inner boundary is located some- 
what closer to the black hole than the other progenitor models. However, as we shall see, 
if the choice of the terminal Lorentz factor is the same, the overall profiles of jet dynamics 
do not depend on whether the jet injection is kinetic dominant or thermal dominant (we 
demonstrate Tj = 5,50 cases in WRLo5 and WRLo50 models). In order to study the de- 
pendence on the accretion-to-jet conversion efficiency, we vary it as r/ = 5 x 10~^, 2 x lO"'^, 
and 10~^ (models such as WRef..., IpopSef..., mpopSef...), while other parameters are iden- 
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tical to the reference model. We also study the dependence on the half opening angle with 
= 3°, 6°, 18°, 36°, and 45°. The study for the injection timing is done only in WR progen- 
itor as WRM3 and WRM6 models. As we have already mentioned, we check the dependence 
on the location of the inner boundary (models such as WRin..., IpopSin..., mpopSin...) and 
we also conduct the resolution checks (models such as WRreso..., IpopSreso..., mpopSreso...). 

According to these studies, we find favorable conditions for creating GRBs. In consid- 
eration of these results, we also perform numerical simulations with representative param- 
eters (models such as WRrepr, IpopSrepr, mpopSrepr). As shown in the following section, 
representative models succeed a powerful relativistic jet breakout and they are the most 
guaranteed candidates to create GRBs in our models. 

3. Result and Analysis 

In this section, we describe the numerical results obtained from our hydrodynamic 
simulations and analyze them in detail. We first explain overall features of the jet dynamics, 
and then we further analyze the dependence on each parameter and model. We also present 
simple analytic criteria for the possibility of GRB production at the end. 

3.1. Basic Feature 

We summarize numerical results in Table [31 Here, we define the shock breakout as that 
the forward shock wave successfully reaches the stellar surface. It should be noted, however, 
that shock breakout is a minimal requirement for producing GRBs. As we shall see, even if 
the outflow successfully accomplishes the breakout, some models are not suitable for GRBs 
(see the column of ^^Possibility of GRB production" in Table [3]). We assess the possibility 
of GRB production by the diagnostic terminal Lorentz factor T^t = h x T profiles on z-axis 
at the time of the shock breakout. (Note that the diagnostic terminal Lorentz factor is the 
achievable Lorentz factor after the internal energy is converted into the kinetic energy.) If 
there are regions where F^j > 100 is satisfied, we determine that the model has a potential 
to produce GRBs. If this condition is not satisfied, the jet is not successfully injected or 
does not move forward. As a result, the explosion would never become relativistic enough. 
Note that some models are hard to be judged for the possibility of GRB production because 
of several reasons (see the column of ^^Possibility of GRB production" in Table [3] and we 
describe these models as A. See also the discussion in Section r3.6p . The details of these 
models are presented in the following subsections. One of the remarkable results in this study 
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is that many models successfully accomplish the relativistic shock breakout even though the 
progenitor star has an e xtremely large envelo pe such as Pop III stars. These results are the 



numerical verification of ISuwa fc lokal (1201 if ). We will further analyze the dependencies on 
the efficiency 77 and the opening angle 9op in the following section. Indeed, the jet dynamics 
strongly depend on these key parameters. 

At the beginning of the simulations, the infall starts at the innermost regions and 
subsequently a rarefaction wave propagates outwards. The accretion rate increases with time 
and then the injected energy generates the forward shock wave around the polar regions. It 
should be noted, however, that the forward shock wave does not propagate outward into the 
active computational regions for a while because of the interruption by the infalling matter. 
If the injected energy is not enough to push the infalling stellar mantle aside, the mantle 
advects inwards and eventually it is swallowed into a black hole. As a result, the total amount 
of explosion energy becomes less than the injected energy. Furthermore, the binding energy 
by a black hole (and also the progenitor star itself) further reduces the explosion energy. We 
also calculate the diagnostic energy E^g at the time of the shock breakout in each model and 
show them in 5th rows at Table [31 We define the diagnostic energy E^g as the integral of 
eic (local energy density), which is the sum of the internal, kinetic and gravitational energy 
density, over the regions with positive eic and . Note that Edg is not the isotropic energy 
but coUimation-corrected true energy. In the weak gravitational field limit, we can define eic 
as; 

eic = - poTc^ + Po^, (4) 

where T** and ip denote the time-time component of energy momentum tensor and gravi- 
tational potential, respectively. Note that, for this definition of diagnostic energy density, 
gravity is treated as an external force. This estimation is valid only when the central core 
object takes the major role for gravity. If the outer envelope plays an important role for 
gravitational energy, we need to take into account the self-gravity. So the term of gravita- 
tional energy should be contributed as a form of not poip but |po^ in Eq. (jl]). However, 
since we use this value for diagnosing whether the outfiow is relativistic or not in this study, 
the severe definition of this quantity is not necessary. Thus, we use Eq. (jll) in the present 
paper. As you can see in Table [3l we find that the diagnostic energy for all models is less 
than about a third of the injected energy. Thus, it implies that the central engine has to 
produce a larger amount of energy than that observed as GRBs. 

The total amount of explosion energy would be larger than the diagnostic energy E^g 
at the breakout since the central engine is able to keep operations after the shock breakout 
due to the continuing accretion. Indeed, as we shall see in the following subsections, some 
models are active at the time of the breakout. Thus, although E^g may be lower than the 
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typical GRB energy in some models, the outflows would gain further energy from the central 
engine to create GRBs after the breakout. 

When the injected energy is successfully launched from the inner boundary, some por- 
tions of matter bounce back and move outwards. However, we can not see a clear coUimated 
outflow at the beginning of the simulations even if we inject kinetic dominant outflows in 
a small opening angle. This is due to the fact that the injected energy is thermalized by 
the strong reverse shock wave in the vicinity of the inner boundary. As a result, the hot 
matter expands and creates a quasi-spherical forward shock wave. Nevertheless, we find 
that the forward shock wave is not strong enough to cease the infall of matter around the 
equatorial region and the matter continues to accrete through the inner boundary (see e.g. 
Fig E] and E]) • The jet structure emerges when the injected energy becomes enough to push 
away the reverse shock wave outside the computational region. Subsequent evolution is sim- 
ilar to the previous studies of jet propagation. The hot cocoons cause recoUimation shock 
waves, and strong backflows also appear in the flows and they sometimes pinch the jet. The 
Kelvin-Helmholtz instability also works to create rich internal structures. The jet starts to 
propagate and the forward shock wave eventually breaks out of the stellar surface. 



3.2. Dependence on the Accretion-to-Jet Conversion Efficiency t] 

As shown in Table |3l we find that the forward shock wave can break out if the accretion- 
to-jet conversion efficiency is 77 > 10~^ (see more details in the subsection 13. 7p . Figure H] 
displays the time evolution of the radius for the forward shock wave on the z-axis for models 
with different efficiencies. As expected, the higher conversion efficiency generates stronger 
shock wave, which quickly propagates into the stellar envelope. In the case of lower conversion 
efficiency models, however, it is harder to inject the jet. Even if the jet is successfully injected, 
it takes a longer time for the forward shock wave to reach the stellar surface than in the high 
efficiency models. 

Figure O shows the time evolution of the accretion rate and luminosity in these models 
(note that we also display the results for the failed models without breakouts). At the 
initial phase, the higher efficiency models have slightly smaller accretion rates than the 
lower efficiency models because the strong outflows interrupt the infall of matter. Still, 
the high efficiency models have higher jet luminosity than the low efficiency models. As 
a result, the jet successfully gains energy from the central engine and accomplishes the 
shock breakout. It is also interesting to note that low efficiency models often show rapid 
time variability of the accretion rate (see e.g. WRef2e-4 model in the upper left panel in 
Figure E]). This is attributed to the fallback of the shocked envelope, which have rich internal 
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structures. Although the envelope initially expands due to the deposited energy, they do 
not have enough energy to keep moving outwards and they eventually fall back through the 
inner boundary, leading to the central engine activity. 

Figure |6] shows the map of regions with positive eic and at the time of the shock 
breakout. We can see that the shape of the yellow region, which contains both positive 
and (see Eq. @), does not depend on the accretion-to-jet conversion efficiency t] so 
much. Irrespective of rj, the maximum transverse radii of the yellow region from the z-axis are 
~ 5 X lO^cm, ~ 2 X lO^^cm and ~ 1.5 x lO^^cm for WR, lpop3 and mpop3 models, respectively. 
It may imply that the outflow structure does not mainly depend on the efficiency t], as long 
as the shock breakout occurs (see also the next subsection). From this figure, we can also 
confirm that even the large efficiency jet does not expel all portions of the stellar mantle, 
and matter can continue to accrete onto the black hole. On the other hand, the inner parts 
of the envelope profiles depend on the efficiency rj. As expected, a larger amount of matter 
is captured for the lower conversion efficiency. In addition, as discussed above, the fallback 
of matter causes the late time variability of the central engine. 

In Figure [TJ we show F^^ profiles along z-axis at the time of the shock breakout for 
each model. The model has the potential to produce a GRB if Tdt > 100 is satisfied in 
some region. Therefore we use this condition to assess whether the model can create GRB 
or not. We note that the reality is more complex as argued in the following. In Figure [TJ 
the outer region tends to contain lower values of T^t, even less than ~ 10. The small F^^ is 
caused by the baryon pollution in the jet, most likely because of the lack of the numerical 
resolution (see in Section I3l6|) . However, even if the baryon pollution were the real physical 
phenomenon, these outflows would create GRBs for the high efficiency models (77 = 10~^) 
because the central engine keeps operating after the shock breakout for these models (see 
discussions in Section IXB]) . The inner fast moving jets would eventually catch up the outer 
slowly moving ejecta. Due to this energy input from the inner jets, it is quite likely that the 
actual terminal Lorentz factor of the outer outflows is larger than the values of the current 
estimation. Thus we expect relativistic explosions in these models. However, it is less likely 
that the low efficiency models such as lpop3ef-4, lpop3ef2-4 and mpop3ef2-4 can successfully 
gain energy from later jets. Even if the central engine keeps operating for a long time, they 
may not be able to accelerate outflows because the amount of the outer slow matter is large 
(the proflle of Tdt deviates from the high efficiency models) in these models. As a result, 
the energy of jet is dissipated in the vicinity of the injection region without transferring the 
energy into the outer parts of the ejecta. According to this consideration, we regard that 
these low efficiency models would end up with non-relativistic explosions and may not create 
GRBs (thus, we mark triangles (A) for these models in the column ^^possibility of GRB" of 
Table 131). We note that if the reverse shock is stalled near the injection region at the time 
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of the shock breakout, there is no hope for this reverse shock to go ahead afterward because 
the accretion rate is decreasing. 

3.3. Dependence on Opening Angle 6'op 

Figure |H] shows the time evolution of the forward shock wave on the z-axis for models 
with different opening angles of the jet at the injection site. Generally, the forward shock 
wave propagates fast and succeeds in breaking out if the opening angle is small. It is mainly 
attributed to two reasons: One reason is that the isotropic flux is increased by the small 
cross-sectional area of the injected jet, and the other is the increment of the accretion rate. 
A wide opening angle tends to interrupt the accretions and reduce the jet luminosity (see 
Figure |9]) . 

Note that the overall dynamics for very wide opening angle jets (e.g. Qop = 36 and 45°) 
are more complex than the collimated jets. As shown in Figure IH1 the jet breakout time 
for 6op = 45° models are (slightly) earlier than 6op = 36° models in all types of progenitors. 
It is because the fall back process plays an important role for wider jets by enhancing the 
accretion rate. Indeed, for WRop45, the accretion rate reaches M ~ 1Mq/s at t = 45 s, 
and leads to the strong outflows from the inner boundary and finally to the shock breakout. 
Thus, the late time activity of the central engine is possible for a wide opening angle due to 
the strong fallback accretions, even though the mean accretion rate decreases with time in 
the late phase. (But these wide opening angle jets are not suitable to produce the relativistic 
outflows, see below). 

In Figure [101 we display the time evolution of the forward shock wave for different angle 
radial rays from the symmetry axis. For reference models (left panels in this figure), the 
forward shock velocity decreases with increasing angle. Particularly, at the time of the shock 
breakout, we find that the forward shock waves along the off- axis radial ray are still deep 
inside the star, which means that the outflow is well collimated. On the other hand, for 
6op = 36° models, almost all models experience the quasi-spherical evolutions. Although the 
shock wave on z-axis is slightly faster than the off-axis shock waves, it is much more spherical 
than the reference model at the breakout. As a result, it is expected that the wide opening 
angle models never create GRBs in view of the spherical morphology of the outflows. 

Figure [TT] is the same as Figure El but for models with different opening angles. We do 
not find any clear differences for narrow opening angle jet models {6op = 3° and 6°), so we do 
not display these models in this figure. As shown in this figure, the wider jet tends to expel 
larger outer envelope, and clearly these profiles are different from each other. Also we again 
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see that very wide angle cases are very complicated. Interestingly, some fraction of stellar 
mantle around the equatorial region may never fall back to the black hole (see 9op > 36° 
models). From what has been discussed above, we can conclude that the outflow profile 
is mainly determined by the opening angle of the jet 6op rather than the accretion-to-jet 
conversion efficiency r] (See also Figure |6] and discussions in subsection 13. 2p . 

The diagnostic terminal Lorentz factor profiles along the z-axis for different opening 
angle models are shown in Figure [T^l For small opening angle models {6op < 9°), the 
diagnostic terminal Lorentz factor is not different very much (slightly larger value for wider 
opening angle). However, the profile of very wide models are completely different from narrow 
jet cases and the terminal Lorentz factor for the wide angle models are quite low. Thus, they 
are no longer capable of producing GRBs. Note that, although the jet of mpop3op36 model 
is successfully injected around the inner boundary (see in Figure [12]) . this outflow may not 
create GRBs since the outer ejecta is completely non-relativistic Tat ~ 1 and the outflow 
configuration is not coUimated. As a result, a small opening angle 6op < 20° is necessary for 
the relativistic shock breakout and GRBs with the case rj ~ 10~^. Note that if the efficiency 
becomes lower than t] = 10~^, the opening angle needs to be smaller than the current value, 
because it is harder for the low efficiency jet to push aside the infall matter than the high 
efficiency jet (see §3.71 for the analytic estimation of the breakout criteria). 



3.4. Dependence on Injection Lorentz factor Tj and Injection Timing tiate 

In this section, we discuss how the injection Lorentz factor and the jet injection timing 
affect the evolution of the jet. In Figure[T31 we display time evolutions of some key quantities, 
such as the forward shock wave on the z-axis, the mass accretion rate and the jet luminosity. 
As we can see, the injection Lorentz factor does not change the qualitative feature of the jet 
evolutions for all quantities as long as the final coasting Lorentz factor is the same (See left 
panels in Figure [13]) . 

On the other hand, the jet dynamics depend on the timing of the jet injection (See right 
panels in Figure [T3]) . We can see that the late injection leads to a shghtly faster evolution 
than the early injection, and also that the time evolution of the accretion rate and the 
jet luminosity are quite different in the early phase. This is because the accretion rate is 
already high at the time of the jet injection for the late injection. As a result, the strong jet 
is injected and evolves faster than the early injection case. It should be noted that although 
the late jet injection slightly interrupts the infall of the material (see middle and right panel 
of Figure [T3]) . it does not cease the accretion of all the infalling material. We also speculate 
that a wide opening angle would suppress the infall of matter more strongly, so that the 
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collimated outflows are preferred for GRBs even in the late injection cases. 

It is also interesting to investigate the diagnostic terminal Lorentz factor r,;;^ profile for 
the different timings of jet injection as displayed in Figure [T^ As we can see in this figure, the 
late jet injection tends to have higher diagnostic terminal Lorentz factor than the early jet 
injection. This is also because the jet luminosity for the later injection case is stronger than 
the earlier jet (see bottom panel of Figure [T3|) . The strong forward shock waves propagate 
outwards and the outer envelope obtains large energy from them. Besides, we also find that 
the late jet injection tends to have a smaller amount of baryon mass in the jet since a portion 
of envelope matter has already been swallowed in the black hole. As a result, the outflow 
easily achieve a relativistic velocity after the breakout (see Table [Hand Section |3T6] for more 
details). We again must caution that the artificial baryon pollution affects the distribution of 
diagnostic terminal Lorentz factor and the amount of baryon mass. The detailed discussions 
of these issues are described in Section 13.61 Note also that too late an injection may not 
be suitable for GRBs, since the accretion rate become s very low and a centrifugal bounce 



takes place at the later time (see iNagakura et al.ll201ll ). Although we do not know exactly 



the starting time of the central engine, we speculate that the most plausible starting timing 
for the central engine is the time when the accretion disk is formed around a black hole. 
It clearly depends on the angular momentum distribution of the progenitor star, so we will 
investigate this dep endence by performin g the rotational collapse of progenitor stars in the 



forthcoming paper ( iNagakura et al.ll2012l ) 



3.5. Representative Models 

Based on the above results, we construct representative models, which are expected to 
create GRBs in each progenitor. As we have seen, the favorable conditions for the relativistic 
shock breakout are the high conversion efficiency, small opening angle and late jet injection. 
We summarize these parameters of representative models in Table |2l We denote these 
representative models as WRrepr for the Wolf-Rayet progenitor, IpopSrepr for the light 
Pop III star and mpopSrepr for the massive Pop III star, respectively. The accretion-to-jet 
conversion efficiency is set as 77 = 10~^ which corresponds to nearly the maximum conversion 
efficiency for a Schwarzschild black hole (See Section |2]). The half opening angle is set as 
^op = 9°. We start to inject the jet when the mass accretion rate becomes the largest for 
each progenitor. We also note that we carry out simulations with higher spatial resolutions 
(7-level AMR) in order to suppress the artificial baryon pollutions (See in Section 13.61 for 
more details.). 

We summarize numerical results in Table El As you can see in this table, the forward 
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shock wave propagates faster than the reference modeL We also find that both the total 
amount of injection energy and diagnostic energy at the breakout time are larger than 
those of reference models. Figure [T5] shows the Tat distribution along the jet axis at the 
shock breakout for each model. As you can see in this figure, every models produces more 
relativistic breakout than the corresponding reference models. In particular, it is interesting 
to note that mpopSrepr succeeds highly relativistic breakout in spite of the large massive 
envelope. Note that the outer ejecta for WRrepr and IpopSrepr is still mildly relativistic 
{Tiit ~ 30). This is most likely caused by the numerical baryon contamination. However, 
even if the baryon pollution were real, these outflows could accelerate relativistically and 
create GRBs for the same reason discussed in Section [3^2] (see also Section [3^ . 



In this subsection, we give some important cautions for the results presented in this 
paper. Although our simplification of numerical methods does not change the essence of 
our new findings, there are some technical limitations for the current numerical simulations. 
Here, we point out several limitations on the present work with some additional simulation 
results. 



Although our numerical code succeeds in capturing strong shock waves and complex 
turbulence in relativistic outflows, numerical diffusion is inevitably inherent in it. As a 
result, numerical diffusion potentially induces artificial baryon pollution which leads to a 
lower diagnostic terminal Lorentz factor in the jet. As we have already mentioned, since 
r^t is important in determining whether the jet would eventually create a strong relativistic 
outflow or not, we need to know how the numerical resolution affects of Tdt along the jet 
axis. 

Figure dn] shows the comparison of F^^ distribution along the jet axis for models with 
different spatial resolutions. As you can see in this figure, models with higher resolution 
exhibit a larger Tat- This confirms that higher resolution reduces the baryon pollution and 
shows the higher value of Tat- In fact, the amount of baryon mass around the jet axis at 
the breakout time is decreasing with higher resolutions (see Table H]). The baryon mass 
contained in the jet is roughly estimated as 



3.6. Limitation of the current study 



3.6.1. Baryon pollution 




(5) 
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where Pzi^) denotes the density profile along the z-axis. We also estimate the average 
terminal Lorentz factor Tj for the jet as 

where Ej denotes the total injected energy after the breakout. We estimate Ej as 

Ej = r]MboundC^ (7) 

where Mhound denotes the total mass of gravitationally bound matter (see also Eq. H] for the 
definition of bound matter). Here, we assume that all of energy from the central engine is 
successfully transferred to outgoing ejecta. As shown in Table HI the outflow for WRref can 
accelerate ~ 60 although this model is heavily influenced by the baryon pollution. Thus, 
even if the baryon pollution were taking place in real, this model could create the relativistic 
outflow. In addition, we would like to emphasize that the macroscopic jet evolution is not 
sensitive to the numerical resolution (see Figure [T7|) . 

We also point out that the two-dimensional axisymmetric setup affects baryon pollution. 
In our simulations, the baryon pollution near the pole is mainly attributed by the numerical 
diffusion, since the meridian velocity near the pole is almost zero due to the axisymmetric 
property. However, in reality, non-axisymmetric motions and finite meridian velocity may 
happen in the jet region. In this case, there is a possibility that the cocoon or strong back 
flows thrust into jets, then a lot of baryons could mix with the jet matter. Note also that the 
dense "plug" at the head of the jet is an artifact of the axisymmetric property. I n reality, the 



non-a xisymmetric motions of jet would cause dispersion of the dense "plug" (IZhang et al. 



20M). 



3.6.2. Dependence on Inner Boundary Location 

One of the other drawbacks in the present study is the choice of inner boundaries in 
current models. For all simulations, the inner boundary is located well outside the central 
core of the star. As a result, the mass accretion rate for each model would be different 
from the actual accretion rate in the vicinity of black hole. In addition, the impact of 
feedback would also be changed. In order to remove these uncertainties, we demonstrate 
some numerical simulations for different position of the inner boundary. 

Figure [18] shows the dependence on the inner boundaries for the evolution of the forward 
shock wave on the z-axis and the jet luminosities. As shown in these panels, models with 
a smaller inner boundary imply that the forward shock wave propagates faster. Since the 
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density is higher at the inner radius, the accretion rate consequently becomes larger so that 
the jet luminosity also becomes larger. In these panels, we find that the model lpop3in5e9 
has a peak accretion rate of about 10 times that of model IpopSref, while WR and mpop3 
simulations exhibit only the factor of 2 or 3 amplification, when the inner boundary is smaller 
by a factor of 2. We also note that, for the mpop3 models, the initial enclosed mass at the 
inner boundary is 414 Mq and it corresponds to nearly half the total mass of the star. Thus, 
the location of the inner boundary substantially affects the jet dynamics especially for lpop3 
and mpop3. In order to obtain more quantitative arguments for the actual jet luminosity 
and breakout time, we need the simulations with a smaller inner boundary, which are beyond 
the scope of this paper. 

It should be noted, however, that the effect of location of the inner boundary is quite 
systematic: the forward shock wave becomes fast and the luminosity becomes large if we put 
the inner boundary on a small radius. This is probably because the density is initially higher 
near the boundary. Consequently the total mass accretion rate increases and therefore the 
jet luminosity becomes large with a strong forward shock wave. According to these results, 
it is a robust claim that the shock breakout is possible if the accretion-to-jet conversion 
efficiency satisfies the condition rj > 10""^ (see more details in the subsection 13. 7p . 



3.6.3. The time delay between infall and injection 

In the present study, we assume that the jet outfiow is injected immediately after the 
mass infiows through the inner boundary. However, it takes a certain amount of time for 
the fiuid to reach the vicinity of the black hole in reality. In addition, the outfiow should 
travel up to the inner boundary. The time difference can be roughly estimated based on the 
assumption that the matter is free fall to a black hole. At the beginning of our simulation, 
it takes ~ 1 s for WR model from the inner boundary to the center of the core, while ~ 20 s 
and ~ 150 s for lpop3 and mpop3, respectively. Our immediate energy injection would be 
different from reality and potentially affect our findings. Here, we investigate the infiuence 
of the time delay between infall and jet injection. 

We perform numerical simulations for some models with taking into account the delay 
time of the jet injection. In this study, we estimate the delay time as free-fall time from the 
location of the inner boundary to the Schwartzschild radius. Note that the Schwartzshild 
radius is estimated by the enclosed mass at the inner boundary. We ignore the jet propagation 
time from the vicinity of the black hole to the location of the inner boundary. This is 
because the time scale of jet propagation is much less than the free-fall time. The jet 
injection parameters, such as the conversion efficiency and opening angles, are the same 
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as the reference models. Hereafter, we denote these models as WRdelay, IpopSdelay and 
mpopSdelay for WR, lpop3 and mpop3 progenitors, respectively. 

In Table 13], we show the summary of results for these models. In comparison with 
reference models, the overall dynamics are qualitatively similar to those of reference models. 
In addition, the variation of the breakout time, the total injection energy, and the diagnostic 
energy remain within 10% from those of reference models. Thus we confirm that the time 
delay between infall and injection is a minor effect. This seems to be attributed to the fact 
that the accretion rate is almost constant except for at the very early phase of collapse. As 
a result, the delayed time injected jet goes through comparable ram pressure to the non- 
delayed injection case. Therefore the forward shock evolution is similar to the reference 
models. 



3.6.4- Long term simulations 

In the current study, we investigate the jet propagations inside of the stellar mantle and 
discuss the possibility of GRBs. However, in order to know how these outflows accelerate 
in the ISM, it is necessary to carry out simulations with longer duration and a large spatial 
region. Although our simulations have a limited duration and spatial extent, we perform 
extended numerical simulations for the WR reference model for the purpose of qualitative 
understanding of the outflow properties. We extend the computational boundary to r ~ 
5Rstar and carry out the simulation until the forward shock wave reaches the outer boundary. 
We cover 1500 uniform radial grids in this extended spatial region, i.e. the total number 
of our radial grids are 2000 (500 + 1500). The grid width is the same as the outermost 
radial width of previous calculations. The outer boundary in this calculation is located at 
i^out = 2.34 X lO^^cm. In this simulation, the 3-level AMR technique is employed. Thus the 
radial resolution in the extended region corresponds to Ar = 4.27 x lO^cm. The meridian 
grid width is the same as previous calculations. 

Figure Uni shows some snapshots for diagnostic terminal Lorentz factor distribution along 
jet axis. As we can see in this figure, the outflow propagates into this ISM and the inner ejecta 
with high Tdt gradually progresses outward with time. This is attributed to the continuous 
energy injection of the central engine. According to these results, even if the baryon pollution 
were the real, the outflows would continuously propagate into the ISM because it is being 
pushed by high Lorentz factor ejecta in the back. Note also that, since this simulation is 
contaminated by the artificial baryon pollution, the outflow is more relativistic in reality 
than this simulation. The encouraging tendencies increase the possibility of GRBs being 
produced in a wider range of scenarios, but future Ion-duration simulations are needed. 
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In summary of this subsection, in order to judge whether the outflows become relativistic 
or not in reahty, we need to simulate them with (1) high resolution in order not to affect 
the baryon pollution, (2) smaller inner boundary, and (3) long duration and large spatial 
range. Here, we emphasize that the results of current study give conservative claims for 
the possibility of relativistic jet breakout, which is the minimum requirement for producing 
GRBs. In fact, we show that the accretion-powered jet succeeds to break out relativistically 
from massive Pop III progenitors if the adequate conditions are satisfied. 



3.7. Comparison with Analytic Estimate 



In this subsec tion, we compare ou r numerical results with the analytic estimate of the 
shock breakout in ISuwa fc loka] (120 111 ) and Appendix IA.2[ Since setups are rather different 
between the numerical study and the analytic one, we concentrate on the parameter depen- 
dence of the successful relativistic shock breakout. They derived the shock breakout time (tb; 
Eq. 12 in their paper) and the free-fall timescale of the envelope materials, i.e., the duration 
of the accretion-powered jet {tg; Eq. 15). When tb is shorter than tg, the jet successfully 
arrives at the stellar surface and the relativistic shock breakout takes place. Employing Eqs. 
(lAlOp and (lAlip in Appendix IA.2t we can derive the criteria for the successful relativistic 
shock breakout as 

6 



— ^ A 



V 
10- 



^ > 1, 



where n 2.6 for mpop3 and WR and ^2.1 for lpop3) is an index for the density profile 
of the stellar envelope, 

pWocf^-lV, (9) 



with the stellar radius i?^, (IMatzner Sz McKedll999l ). The pref actor A depends on the stellar 
radius, the core mass a nd the envelope ma ss, and is found to be typically around A ~ 3 
for mpop3 according to ISuwa fc lokal (120111 ) and Appendix 1A.21 In this paper, we find that 
A ~ 2 is suitable for explaining the numeri cal results of mpop3. The reason of the difference 
in A is that some of the assumptions in ISuwa fc lokal (12011! ) are violat ed, that is, they 
assumed a conical jet propagation (i.e., 9op is independent of the radius; see lloka et al.ll201ll ) 
and neglected the negative feedback from the cocoon on the accretion rate. The breakout 
condition is also not exactly the same between the analytic and numerical calculations. 
Nevertheless it is encouraging that both results coincide within a numerical factor. 

In Figure [20] we present r] — 6op diagrams to compare our numerical results with the 
analytic criteria in Eq. (black solid line). The circles, triangles and crosses in this figure 
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correspond to the column of ^^Possibility of GRB production" in Table |3l We use A = 2 
for all progenitor rnodels. As we can see from Figure |20l the analytical criteria in Eq. ([8]) 
and ISuwa fc lokal (120 111 ) are useful to forecast the penetrability of the relativistic jet and 
the possible GRB production in massive star models. We can predict the parameter space 
where GRBs could occur in — 9op plane once we calibrate A with numerical simulations. In 
addition, we find that the analytical critical curve can be approximated by a much simpler 
form for all current progenitors (WR, mpop3, lpop3) as 



^> 10- 



op 



(10) 



Note that a red supergiant (RSG) satisfies the breakout criteria in Eq. (|8]) with A ~ 6.8. 
Actually a GRB jet might be associated with a RSG, which has not been observed so far. 
A GRB from a RSG would be very dim and long because the breakout time is long due 
to the expanded envelope and the accretion rate is low at the late time. Such dim GRBs 
might even dominate as the sensitivity is improved. Alternatively, another condition could 



(Matzner 


2003: 


Suwa & loka 


2011) 



In some parameter space, a 
RSG does not satisfy this second criterion be cause of its shallow envelope and the outflow 



becomes almost spherical ( ISuwa fc lokal 1201 ll ). Second, the central engine of a RSG might 



not rotate fast enough for the jet production. The detailed discussion will be presented in 
the forthcoming paper. 



4. Conclusions 

In this paper, we investigate the propagation of the accretion-powered jets in Pop III 
and present-day stellar progenitors (Wolf-Rayet stars). We perform two-dimensional rela- 
tivistic hydrodynamic simulations taking into account both the envelope collapse and the jet 
propagations for the first time. The main findings in this paper are summarized as follows. 

1. Although some portion of matter ceases to fall and is pushed outward by the energy 
injection from the polar region, a large amount of matter continues accreting into the 
black hole, so that the central engine can continue to inject outflows. 

2. If the central engine satisfies a certain condition (see below), the jet can successfully 
propagate and create a relativistic breakout from various types of progenitors for GRBs. 
In particular, we show that Pop II I stars could be progenitor candidates for GRBs, as 



pointed out by ISuwa fc lokal (120111 ) . We numerically verify that the central engine can 
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last very long > 100 s for light Pop III stars and > 1000 s for massive Pop III stars 
because of the accretion of the huge envelope. 

3. The jet can produce a relativistic breakout if the accretion-to-jet conversion efficiency 
satisfies 

^>io-(^)', (11) 

as derived in Eqs. ([8]) and (ITO!) . Otherwise the injection energy results in non-relativistic 
explosion or advection into a black hole. 

4. We find that the timing of jet injection does not affect our results significantly. We 
also observe that the diagnostic terminal Lorentz factor tends to be higher for later 
jets than for early jet injection because the rarefaction wave reduces the central density 
before the strong jet is injected. 



It should be noted that our simulations are affected by numerical diffusion and the 
location of inner boundary as analyzed in Section 13.61 Especially for Pop III models, we put 
the inner boundary far outside the central core. We intend to address these issues in future 
publication. 

We would like to point out that the injected energy from the central engine may con- 
tribute to the explosion energy of the supernova component, if the hot cocoon mainly con- 
tributes to the supernova, i.e., an almost spherical s h ock w ave accompanying a large amount 



of nickel production (See also iRamirez-Ruiz et al.l (120021 ) for discussions of excess energy 



accumulated in the cocoon). If the strong reverse shock wave is still stagnated around the 
root of the jet at the time of the shock breakout, the kinetic energy of the jet continues 
to be converted into thermal energy deep inside of the stellar mantle. Although the cocoon 
pressure would decrease after the shock breakout and make it easy for the reverse shock wave 
to propagate, the jet luminosity is also reduced as the accretion rate decreases, so that the 
reverse shock wave may stay deep inside the star. If this is the case, the injected energy may 
work to expel the stellar mantle rather than contribute to the GRB component. However, 
at the present time, it is unclear how large an amount of nickel is created in the hot cocoon, 
so we do not know whether the cocoon creates the supernova or not. We will address these 
issues in a forthcoming paper. 

It is also interesting to note that if the stellar envelope consists of multi-layerd config- 
urations (i.e. an onion like structure), the accretion rate also changes across the different 
layer. The change of the accretion rate may cause time variability of the central engine. The 
variability timescale depends on the initial location of the shell boundaries and roughly cor- 
responds to the free-fall timescale at this region. It should be noted that the stellar rotation 
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may play an important role on this matter for the realistic situation since the specific angular 
momentum is also discontinuous at the shell interface, and it may cause the accretion rate 
to fiuctuate. The angular momentum is generally an increasing function of radius in each 
shell, but it discontinuously drops at the shell boundary. As a result, the centrifugal force is 
reduced in strength there and it may potentially increase the accretion rate. 

It is interesting to note that some models in the present study may explode as choked 
GRBs. Although they can n ot produce GRBs, they are expected to produce TeV neutrinos 
( iMeszaros fc WaxmanI l200ll ). These observational signals provide important information 
on jet propagation in the stellar mantle. We would like to investigate the observational 
consequence of the different jet injection parameters in the forthcoming paper. 

Finally, we would like to emphasize that the rotation profile of the stellar progen- 
itor also affects the GRB production. This is because the angular momentum profile 
may determine not only the starting time of the central engine but also the disk cond i- 
tions ( iLee &: Ramirez-Ruizll2006l : IZalamea fc Beloborodovll2009l : iLopez-Camara et al.ll2009l ). 
The dependence on an gular momentum profile of progenitor is currently being investigated 
( iNagakura et al.ll2012l ). Additionally, there is no guarantee that the jet luminosity is propor- 
tional to the accretion rate as assumed in this paper. Actually, if the neutrino mechanism 
plays a key role for producing the jet, the jet luminosity also depends on the mass of the 
black hole. In addition, the relation between the mass accretion rate and jet luminosity 
in neutrino mechanism would be different from those used in the present paper (see e.g. 
Zalamea fc Beloborodovll201ll ). These areas present the opportunity for future inquiries and 
are currently being investigated. 
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A. Appendix 



A.l. Analytic stellar model 



Using Equation (9) of ISuwa fc lokal (1201 ll ). we can calculate the stellar radius as 



= 10^3 cm 




Mr 



IQio cm J V4OOM0 



1 

3 — 71 



V SOOMp 



1 

3 — 71 



(Al) 



where we calibrate the overall factor to repro duce the mpop3 ra dius with n = 2.6. This 
equation is a general form of Equation (10) of ISuwa &: lokal ( 1201 ll ). From Eq. (lAip . we can 
estimate the stellar radii of WR, lpop3, and RSG as 8 x 10^'^ cm, 7 x 10^^ cm, and 2 x 10^'^ cm, 
respectively, which well reproduce actual values within a factor of 2-3. In these estimations, 
we employ following parameters; 



• WR: n = 2.6, = HMq, Menv = 2Mq, and Rc = 10^° cm. 

• lpop3: n = 2.1, = ISM©, Menv = 25M0, and R^ = 10^° cm. 

• RSG: n = 1.5, = AMq, Me„v = 8Mq, and R^ = 10^^ cm. 



A. 2. Stellar dependence of breakout criteria 

In this section, we derive analytic expression for the breakout criteria. Note that we 
employ typical sets of parameters to normalize physical quantities (e.g., rj, Bop). In order to 
apply to the individual progenitor, one should insert the values in Appendix lA.li 

Here, we approximate the density profile of the envelope as 

Pir)=Pi(-]\ (A2) 



in which the term of unity in Eq. (j9]) is neglected for simplicity. This approximation is valid 
for r <C -R*. The free-fall timescale of a mass shell at the radius r and the enclosed mass Mr 
is given by 



Lj. 



where G is the gravitational constant. The mass accretion rate is given by 

dMr/dr ^ Stt 
dtff/dr ^ 3 tff 
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where Rq is the arbitrary characteristic radius and po = pi(-R*/-Ro)"'- Here we neglect the 
derivative of Mr with respect to r for dtff/dr. Using Eq. flA4p . the jet luminosity emitted 
from the central object is written as 

L, = ,Mc^ = 7.5 X 10- erg s^ (^) ( ^V' f^— ( ( /° (A5) 

^ ' ^ vio-v\^isy Vio^icm/ V^/ vio-^gcm-3y^ ^ 

where tff corresponds to the jet breako ut time at -Rn. On the other hand, the necessary 
luminosity for outgoing jet propagation (ISuwa &: Iokall201ll ) is given by 

lJ-^ = 1.3 X 10^« erg s"^ ( -4^) ' ( ,1 V f ^ V' . (A6) 

Equating Eqs. f lASP and ( ]A6p . we get 



By deleting r from Eq. (jAT]) using Eq. (jM]) (i.e., r ~ 10^^ cm(M^/25OM0)i/3(%/17O s)^/^), 
we can estimate the jet breakout time (i.e., the jet arrival time at i?o) as 

6 3(4-n) 3-Ti 

t.(r = i?o) ~ 170 s (^j (^-J [j^TT^) [-^^) ■ (A8) 
On the other hand, the free-fall time in Eq. (1A3P is 

Here, we choose Rq ~ 0.3-R*, beyond which the density profile is not power law but decreases 
rapidly due to the term of unity in Eq. IQ. The jet head accelerates due to decreasing ram 
pressure of the ambient matter beyond this radius. In addition, the mass accretion rate by 
the mass shell at r > 0.3i?,, decreases rapidly. Therefore, we consider the critical condition 
for the successful jet breakout using Rq ^ O.SR^, in the following. 



By combining Eqs. f lA8|) and f lA9|) . we obtain the criterion for the successful jet break 
out, which is that the free-fall time scale should be longer than the breakout timescale, as 

According to the notation of Eq. ([H]), we define A as 



1011 cmy VSSOMey 
As for Mq,^r , we replace this term with Mc + 0.4Mcnv to estimate the specific values (see 



Suwa fc Iokall201ll ). Using values in Appendix IA.lt we find that A ~ 1.5, 2.6, 2.7, and 6.8 



for WR, lpop3, mpop3, and RSG, respectively. 
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Table 1. Summ ary of stellar models: (1) Wolf-Rayet star (16TI in (IWoosley fc Hegei 
star f lWoosley et al.ll2002[ ) (lpop3), (3) massive Pop III star ( lOhkubo et al. 



2mm, WR), (2) light Pop III 
mm (mpopS). 
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Table 2. Summary of our models 
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Table 2 — Continued 



Model 


Inner boundary 
Ri„ (cm) 


Enclosed mass 
Mi„{M0) 


Efficiency 
V 


Half opening angle 
Sop (°) 


Lorentz factor 


Specific internal energy 
ej 


Retarded injection time 

tlate (s) 


AMR level 


mpopSref 


10" 


414.4 


10-3 


9 


400 


10-2 





3 


mpop3ef5e-4 


IQii 


414.4 


5 X 10-'' 


9 


400 


10-2 





3 


mpop3ef2e-4 


IQii 


414.4 


2 X 10-" 


9 


400 


10-2 





3 


mpop3efle-4 


IQii 


414.4 


10-4 


9 


400 


10-2 





3 


mpopSopS 


10" 


414.4 


10-3 


3 


400 


10-2 





3 


mpop3op6 


10" 


414.4 


10-3 


6 


400 


10-2 





3 


mpopSoplS 


10" 


414.4 


10-3 


18 


400 


10-2 





3 


mpop3op36 


10" 


414.4 


10-3 


36 


400 


10-2 





3 


mpop3op45 


10" 


414.4 


10-3 


45 


400 


10-2 





3 


mpop3in5el0 


5 X 10^° 


385.1 


10-3 


9 


400 


10-2 





3 


mpop3reso7 


10" 


414.4 


10-3 


9 


400 


10-2 





7 


mpop3repr 


10" 


431.2 


10-2 


9 


400 


10-2 


246.6 


7 


mpop3delay 


10" 


414.4 


10-3 


9 


400 


10-2 





3 



Table 3. Summary of our models 
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Table 3 — Continued 



Model 
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^ O-Success of shock breakout, XiFailure of shock breakout. 

^ O-There is a possibility of GRBs, AiTouchy to judge the possibility of GRBs, X:There is no possibilities of 
GRBs. 

Time at shock breakout. This time is measured from the beginning of central engine. 



Table 4. The Baryon Mass and Final Lorentz Factor for Jets 



Model 


Baryon Mass^" {A/0) 


Jet Energy'' (Ej) 


Average terminal Lorentz factor'^ (T/) 


WRref 
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9.17 X 1052 


2035 



^ The baryon mass contained in the jet at the breakout. Sec Eq. (O in Section 13.61 for its 
definition. 

^ The total energy of injected jet after the breakout. 

The average terminal Lorentz factor of outflowrs estimated by Eq. JT)!- 
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The density profiles of the progenitor models used in this paper (see Tabled]). 
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Fig. 2. — The radial density profile for the late time injection models. The red line is the 
original stellar profile of WR model while the green (blue) line indicates the density profile 
at the time when the enclosed mass at the inner boundary becomes M = 3Mq (M = QMq) 
as a result of the spherical collapse. The rarefaction waves can be clearly seen in this figure. 
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Fig. 3. — The time evolutions of inner core mass (black hole mass) for each model. 
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Fig. 4. — Time evolutions of forward shock waves on z-axis for different accretion-to-jet 
conversion efficiency t] models. Left; WR models. Middle; lpop3 models. Right; mpop3 
models. Note that at the beginning of simulations, the artificial oscillations arise due to the 
difficulty of identification of shock position. However, they do not affect our main results. 
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Fig. 5. — The time evolution of accretion rate (left) and luminosity (right) for different 
accretion-to-jet conversion efficiency rj models. The upper, middle and bottom panels show 
WR models, lpop3 models and mpopS models, respectively. 
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Fig. 6. — The regions with positive eic and (Yellow) at the time of the shock breakout 
for different accretion-to-jet conversion efficiency rj models (see section [HH] for the definition 
of eic and f''). The upper, middle and bottom panels show WR models, lpop3 models and 
mpopS models, respectively. 
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Fig. 7. — The diagnostic terminal Lorentz factor profile {Tdt = h x T, i.e., the Lorentz 
factor which the jet can in principle attain) at the time of the shock breakout for different 
accretion-to-jet conversion efficiency t] models. From left to right, WR models, lpop3 models 
and mpop3 models, respectively. 




Fig. 8. — Time evolutions of forward shock waves on z-axis for different opening angle 
9op models. From left to right, we show WR models, lpop3 models and mpop3 models, 
respectively. 
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Fig. 9. — The time evolution of the accretion rate (left) and the luminosity (right) for 
different opening angle 6op models. The upper, middle and bottom panels show WR models, 
lpop3 models and mpop3 models, respectively. 
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Fig. 10. — The forward shock evolution along each radial ray with different angle from the 
axisymmetric axis. We show two different opening angle jet models. Left: reference model 
{9op = 9°), Right: {9op = 36° model). From upper to lower panels, we show WR, lpop3 and 
mpopS models, respectively. 
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Fig. 12. — Same as the Figure [7] but for different opening angles. 
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Fig. 13. — The dependence on the injection Lorentz factor and the timing of the jet injection. 
Left: different injection Lorentz factor, Right: different injection timing. Upper to lower; 
the time evolutions for the shock evolution along the z-axis, the mass accretion rate and the 
jet luminosity, respectively. 
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Fig. 15. — Same as Figure [71 but for representative models. 
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Fig. 16. — Same as Figure [TJ but with different spatial resolutions among WR models 
(WRref, WRresoS and WRreso7). 
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Fig. 17. — The dependence on the resolutions of our simulations. Left: the forward shock 
evolution along the z-axis, Right; the time evolutions of the luminosity. 
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Fig. 18. — The dependence on the location of the inner boundary. Left: the forward shock 
evolution along the z-axis, Right; the time evolutions of luminosity. 
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Fig. 19. — Same as Figure [TJ but long term simulations for WRref. The red line indicates the 
diagnostic terminal Lorentz factor distribution soon after the shock breakout, while the blue 
line indicate the same distribution but t = 14.5 (s) which is the final time of this simulation. 
The snapshot for green line {t = 11.4(s)) corresponds to the middle time between them. 
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Fig. 20. — The score sheet of the shock breakout for WR (top panel), lpop3 (middle panel), 
and mpop3 (bottom panel), respectively, in the accretion-to-jet conversion efficiency rj and 
the opening angle 6op plane. Red circles correspond to the models which have possibilities for 
creating GRBs in our numerical simulations, while blue crosses show the failed cases. Black 
triangles are marginal models (see text for details). The analytical criteria for the shock 
breakout in Eq. are shown by the black solid lines. Below this line the shock wave can 
break out of the stellar surface before the mass accretion ceases. On the other hand, the jet 
stalls inside the massive envelope and the explosion becomes spherical above this analytical 
line. 



